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Adsorption of or^/io-deuterium and j^ara- hydrogen films on a graphite substrate, pre-plated with 
a single atomic layer of krypton, is studied theoretically by means of quantum Monte Carlo sim- 
ulations at low temperature. Our model explicitly includes substrate corrugation. Energetic and 
structural properties of these adsorbed films are computed for a range of hydrogen coverages. Ther- 
modynamically stable adsorbed films are solid, with no clear evidence of any liquid-like phase. 
Quantum exchanges of ort/io-deuterium and j^ara- hydrogen are essentially absent in this system, 
down to zero temperature; consequently, this system displays no superfluidity in this limit. Our 
simulations provide evidence of a stable domain wall fluid at low temperature, consistently with 
recent experimental observations. 

PACS numbers: 02.70.Ss,67.40.Db,67.70.+n,68.43.-h. 



I. INTRODUCTION 

Studies of low temperature adsorption of highly quan- 
tal fluids, such as helium or molecular hydrogen are mo- 
tivated by the fascinating properties that such adsorbed 
quantum films display, often considerably different than 
those of the bulk. Fluids of or^/io-deuterium (0-D2) and 
para- hydrogen (P-H2) molecules are interesting physical 
systems, because these molecules have mass equal to (in 
the case of 0-D2) or half (in the case of p-B.2) the mass 
of a helium atom, rendering zero-point motion signifi- 
cant at low temperature. In addition, each molecule is 
a composite boson, and therefore it is conceivable that, 
at low enough temperature, a hydrogenic fluid might dis- 
play physical behavior similar to that of liquid helium, 
including superfluidity (SF).^ 

However, unlike helium, bulk hydrogen solidifies at low 
temperature (T^ ^ 14 K); this solidification has so far 
prevented the observation of phenomena such as Bose 
Einstein condensation (BEC) and SF which, in the bulk, 
are speculated to occur in the liquid phase below T ^ 

6 K. Solidification occurs because of the depth of the 
attractive well of the potential between two hydrogen 
molecules, which is significantly greater than that be- 
tween two helium atoms. Several attempts have been 
made^i^i^i^ to supercool bulk liquid P-H2 , but the search 
for SF (in the bulk) has so far not met with success. 

Toward the aim of maintaining a stable liquid phase of 
molecular hydrogen down to low enough temperatures for 
SF to be observed, several avenues can be explored, in- 
cluding reduction of dimensionality. The predicted^ the- 
oretical value of the melting temperature for a solid 
film of P-H2 in two dimensions (2D), is slightly less than 

7 K, i.e., approximately half that of three-dimensional 
bulk P-H2, but still significantly higher than that at which 
the system, if it remained a liquid, would undergo Bose 
Condensation and turn superfluid, estimated at ~ 2 K 



in 2D. However, in an actual experimental realization 
of a quasi-2D system, i.e., a film of hydrogen molecules 
adsorbed on a substrate, quantum zero-point motion of 
adsorbed particles in the direction perpendicular to the 
substrate can be significant, as calculations for adsorbed 
helium films on alkali metal substrates have shown. ^ One 
might speculate that such zero-point motion may result 
in an effective screening of the hydrogen intermolecular 
interactions, possibly leading to a reduction of T^, with 
respect to the purely two-dimensional case. 

With this aim in mind, extensive experimental studies 
of hydrogen films adsorbed on a different substrates have 
been carried out, using various techniques ^^i^ii^iiiii^ii^ii^ 
One of the most remarkable findings, is that the melting 
temperature of a solid p-B.2 monolayer^ ^ can indeed 
be significantly less than that of purely 2D P-H2. The 
interesting question is whether by "fine tuning" of the 
substrate, Tm may be lowered to the point where inter- 
esting, collective quantum many-body phenomena could 
become observable, in some thermodynamically stable 
liquid phase. 

In a recent neutron scattering investigation of 0- 
D2 films adsorbed on a krypton pre-plated graphite 
substrate^ evidence of a stable "liquid-like" phase of 0- 
D2 down to T ~ 1.5 K was reported (termed a "domain- 
wall" fluid by the authors). This result is obviously of 
great interest, as it appears that the substrate consid- 
ered in this work may indeed provide an avenue to the 
observation of the long sought superfluid phase of molec- 
ular hydrogen, given the low temperature down to which 
a liquid-like phase may exist. Motivated by this experi- 
ment, we have undertaken a theoretical study of the low 
temperature phase diagram of P-H2 and 0-D2 films ad- 
sorbed on such a substrate. 

In a previous study^ we determined that P-H2, ad- 
sorbed upon a substrate consisting of graphite preplated 
with a commensurate monolayer of krypton forms two 
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t her mo dynamically stable monolayer phases at low T, 
both solid; one is commensurate with the krypton layer, 
the other incommensurate. No evidence was observed of 
a liquid phase at T = for these or any intermediate cov- 
erages, and quantum exchanges of P-H2 were found to be 
greatly suppressed. Quantum zero-point motion in the 
direction perpendicular to the substrate was found to be 
less significant than that observed for P-H2 adsorbed on 
lithium, from which one can conclude that this mecha- 
nism ought not yield a significant reduction in the melt- 
ing temperature. In that study, while the Kr layer was 
modeled explicitly, the underlying graphite substrate was 
assumed smooth. 

In this paper, we again explore this problem, but with 
two main modifications with respect to Ref. [l^, namely 
a) we assume this time a compressed krypton layer pre- 
plating the graphite, i.e., incommensurate with the un- 
derlying substrate lattice structure, and b) we explicitly 
model corrugation of the underlying graphite substrate, 
in some of our simulations. Furthermore, in order to 
allow for a direct comparison with Ref. 15, we study 
adsorption of both P-II2 and 0-D2. Energetic, structural 
and superfiuid properties are investigated using Path In- 
tegral ground state (PIGS) Monte Carlo simulations, as 
well as the continuous-space worm algorithn>i2iii^ for sim- 
ulations at finite temperature. 

The main results of this study are the following: 

1. No evidence is observed of a thermodynamically 
stable liquid phase in 0-D2 or P-H2. In all cases, 
equilibrium phases are solid monolayers, commen- 
surate with the underlying krypton plating, with 
coverage (i.e., 2D density) l9o=0.07253 A-^ for both 
isotopes. 

2. An 0-D2 monolayer composed of commensurate do- 
mains separated by domain walls is found to be 
thermodynamically stable, compressible 

in a relatively extended range of coverages. 

3. Quantum exchanges of p-B.2 or 0-D2 molecules are 
essentially absent in the T ^ limit; consistently, 
no evidence of a finite superfiuid response of either 
isotope is ever observed. 

4. Our results are consistent with the presence of 
a "domain- wall fiuid", as suggested in Ref. [l5|. 
Though no dynamical information can be reli- 
ably extracted using the QMC techniques employed 
here, visual inspection of instantaneous configura- 
tions generated by our sampling procedure suggests 
that such a domain- wall fiuid indeed form; however, 
such a fiuid is classical in character, and does not 
support SF of either hydrogen isotope even at T=0. 

The remainder of this manuscript is organized as fol- 
lows: Sec. nil offers a description of the model used for 



our system of interest, including a discussion of the po- 
tentials and the justifications for the main underlying 
assumptions. Sec. IIIII involves a brief discussion of the 
computational techniques and specific details of imple- 
mentation, in addition to details of calibration and opti- 
mization. The results are presented in Sec. IIV( finally. 
Sec. [Vlis a summary of the findings and our concluding 
remarks. 



II. MODEL 

We consider a system of TV hydrogen molecules, re- 
garded as point particles of nuclear spin zero (i.e., obey- 
ing Bose statistics; this discussion is entirely equivalent 
for both isotopes). Molecules sit above a substrate con- 
sisting of a single atomic layer of krypton, below which is 
a graphite substrate. The Kr layer consists of (L) point- 
like atoms, pinned at fixed positions {k = 1,2, ...,L) 
(an assumption justified by their relatively large mass). 
They are arranged in a perfectly two-dimensional trian- 
gular lattice, with a spacing^^ of 3.99 A and sit at a 
height of 3.46 A over the top layer of C atoms of the 
graphite substrate; this latter distance corresponds to 
the energetic minimum of the most accurate Kr-graphite 
potential available ^^SaSI if graphite is regarded as smooth. 

If a smooth graphite substrate is assumed, the model 
quantum many-body Hamiltonian is therefore as follows: 



i<j 
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Here, m is the mass of a hydrogen molecule, {r^} (with 
z=l,2,...,7V) are the positions of the hydrogen molecules, 
rij = \ri — Tjl; Zi is the height of the ith hydrogen 
molecule above the graphite surface. V is the poten- 
tial describing the interaction between any two hydrogen 
molecules, and U represents the interaction of a hydrogen 
molecule with a Kr atom. Finally, U represents the inter- 
action of a hydrogen molecule with the smooth graphite 
substrate. 

As we are investigating P-H2 and 0-D2, which have 
nearly spherically symmetric electron densities, all pair 
potentials are assumed to depend only on relative dis- 
tances. The interaction V is described by the Silvera- 
Goldman potential^H which provides an accurate de- 
scription of energetic and structural properties of con- 
densed P-H2 and o-B2^^ 

The interaction of a hydrogen molecule and a Kr atom 
is modeled using a standard 6-12 Lennard- Jones (LJ) 
potential; we make use of the Lorentz-Berthelot mixing 
rule^i2^ yielding e = 75.6 K and a = 3.3 A, for our 
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purposes consistent with current state-of-the-art poten- 
tial energy surfaces (for example, numerical agreement 
is within ^ 5% for the well depth) The Kr atomic 
monolayer is uniform, with each Kr atom kept fixed in 
its lattice position for the duration of the simulation. 

We use a simple "3-9" potential to describe the inter- 
action of hydrogen molecules with the smooth graphite 
substrate, i.e.. 



C 



(2) 



where C=7913.24 K and L>=259.39 K are param- 
eters derived from the original P-H2-C Lennard- Jones 
parameters'^ (<j = 3.18 A, e = 32.05 K) and the den- 
sity of carbon atoms in graphite {p = 0.114 A~^). 

The assumption built in ([T]) of a smooth graphite sub- 
strate seems justified, on account of the relatively large 
(~ 6 A) distance between the hydrogen molecules and 
the substrate, due to the presence of the Kr layer. Nev- 
ertheless, corrugation may be important in this problem, 
given that the Kr layer is not registered with the under- 
lying graphite substrate. Therefore, in order to gauge 
the effect of substrate corrugation, we performed a num- 
ber of simulations with an explicitly modeled top layer 
of C atoms, placing a smooth graphite slab underneath. 
For these simulations, the pair- wise interaction between 
a hydrogen molecule and a C atom was taken to be the 
same LJ potential utilized to infer the "3-9" interaction 
of hydrogen molecules with a smooth substrate. 

The results obtained using this composite substrate in- 
dicate no qualitative change in the physical character of 
the system as a result of introducing the smooth plane 
approximation for graphite. Quantitatively, the poten- 
tial energy per particle, for example, was found to be 
^2 K lower when assuming explicit corrugation of the 
graphite substrate, a difference of less than 1%. Hence- 
forth, therefore, we shall confine our discussion to the 
case of a smooth graphite substrate, i.e., model ([1]). 

A slice of the effective potential energy surface for 
molecular hydrogen resting at its average height above 
the substrate is shown in Figure [TJ included are the in- 
teractions with both the Kr monolayer and the graphite 
slab. There are two energetically degenerate sublattices 
of preferential adsorption sites for this geometry - these 
are denoted sublattice 
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and sublattice 
The model ([1]) clearly contains important physical sim- 
plifications, such as the neglect of zero-point motion of 
Kr atoms, as well as the restriction to additive pairwise 
interactions (to the exclusion of, for example, three-body 
terms), all taken to be central, and the use of the highly 
simplified LJ and "3-9" potentials. Nonetheless, it seems 
a reasonable starting point, and even quantitatively we 
expect it to capture the bulk of the physical picture. 




FIG. 1: (color online) Potential energy surface, at the verti- 
cal equilibrium distance of molecular hydrogen, as a function 
of position. The two energetically degenerate sublattices of 
preferential adsorption sites, A and B, are shown. A com- 
mensurate coverage is defined as the occupation of one entire 
sublattice (also termed "| filling"). This is a low resolution 
image; contact the authors for high resolution image. 



III. COMPUTATIONAL METHOD 

Accurate ground state expectation values for quantum 
many-body systems described by a Hamiltonian such as 
^ can be computed by means of Quantum Monte Carlo 
(QMC) simulations. In this work, the method utilized 
is Path Integral ground state (PIGS), which is an exten- 
sion to zero temperature of the standard. Path Integral 
Monte Carlo method*^ PIGS is a projection technique, 
which filters the exact ground state wave function out 
of an initial trial state. It is therefore closely related to 
other ground state projection methods, such as Diffusion 
Monte Carlo (DMC), but has a few distinct advantages 
(for a discussion, see, for instance, Ref. [31I ). Because this 
method is described in a number of publications, it will 
not be reviewed here. Some of the technical details of 
the calculation performed in this work (mainly, the short 
imaginary time propagator) are the same as in Ref. [32I . 

The trial wave function utilized is of the Jastrow type: 



^T(ri,r2,...r^) = (11^"'^'"^ (H 
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X 'nn^ 

i=i 
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with pseudo-potentials w (p-H2-graphite and 0-D2), u {p- 
H2-Kr and o-D2-Kr), and v (P-H2-P-H2 and 0-D2-0-D2) 
chosen as follows: 



w{r) 



a 
73 



u{r) 
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and v{r) 



p_ 



(4) 



The values of the parameters a = 30 A^, 7 = 250 A^ 
and /i = 750 A^ were obtained empirically, by minimizing 
the energy expectation value computed in separate vari- 
ational calculations. Using the trial wavefunction as de- 
fined above, we observe convergence of the ground state 



4 



energy estimates with a projection time 0.250 K~^, using 
a time step r = (z^/2560) K~^, where = 1 for P-H2 and 
= 2 for 0-D2. 

Finite-temperature results were obtained making use 
of a reformulation of standard path integral Monte Carlo 
(PIMC) called the worm algorithm (WA), a scheme 
which allows for efficient computation of thermodynamic 
properties, including winding numbers and off-diagonal 
correlations, for systems of much greater size than that 
accessible to conventional PIMC simulations.—!^ The 
only input parameter to the WA, as for any finite tem- 
perature method, is the microscopic model ([T]). Just as 
in ground state simulations, an optimal value of the 
imaginary time step must be determined, such as accu- 
rate and unbiased estimates of the observables are ob- 
tained; for this system, satisfactory results are yielded 
by the choice r^^ = (i^/320) K~^. 

Calculations for a range of P-H2 and 0-D2 coverages 
were carried out, starting from an initial configuration of 
molecules sitting atop the Kr layer. The simulation cell 
consists of a 12x14 triangular lattice of Kr atoms with 
3.99 A nearest neighbor spacing (resulting in a simula- 
tion box of dimension 47.880 A x 48.3762 A).^ Periodic 
boundary conditions are used in the three directions, but 
the simulation cell is chosen sufficiently large in the z di- 
rection that they have no effect vertically, because of the 
strongly attractive character of the composite substrate 
(i.e., P-H2 and 0-D2 molecules do not evaporate). It is 
important to note that to capture the features observed 
experimentally in Ref. 15 requires that a large number of 
particles be simulated (the simulation cell must be large 
enough to support at least two crystalline domains). 

The systematic errors of our calculation are at- 
tributable to finite projection time (for the T=0 PIGS 
calculations) and the finite time step r (for both PIGS 
and the WA). Based on comparisons of results obtained 
from simulations with different values of projection time 
and/or time step, we estimate our combined systematic 
error on the total energy per P-H2 molecule to be of the 
order of 0.7 K or less (corresponding to less 0.5%). The 
simultaneous use of a ground state and a finite temper- 
ature method allows us to obtain an independent check 
of our calculations in the T ^ limit, where the two 
methods must yield the same results, within statistical 
uncertainties. 



IV. RESULTS 

Physical quantities of primary interest, for both P-H2 
and 0-D2, include the energy per molecule, e(^), and the 
superfluid fraction, pg. PIGS (T=0) was used to compute 
e{0) as a function of the coverage (two-dimensional den- 
sity), 6>, from 19=0.06735 A-^ to (9=0.09498 A-^. For a 
subset of these coverages, the worm algorithm (WA) was 



employed to yield finite temperature estimates for all of 
the above physical quantities, for a set of temperatures 
between T=5 K and T=0.25 K. 
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FIG. 2: (color online) Energy per molecule e{0) (in K) for p- 
H2 (o), shifted by —50 K for clarity, and 0-D2 (•) computed 
by PIGS (T=0), as a function of the coverage (in A~^). 
Inset: e{0) calculated using WA for 0-D2 is shown with T 
for 6'=0.08389 A"^ and V for 6>=0.07253 A"^ values marked 
with □ on the T—0 0-D2 curve. 

The results for e{0) are shown in Figure [21 note that 
e{0) for P-H2 has been shifted by — 50 K for clarity. In the 
case of P-H2, the main feature is an energetic minimum 
at ^=0.07253 A~^, corresponding to the formation of a 
thermodynamically stable solid layer, that is commensu- 
rate with the underlying krypton lattice. This coverage 
is approximately 8% greater than the equilibrium density 
of P-H2 in strictly 2D.^ Within the precision of the calcu- 
lation, e{0) is numerically consistent with the absence of 
other thermodynamically stable coverages, in the range 
explored in this work. 

For 0-D2, we find that there is an energetic minimum at 
19=0.07253 A-2, again corresponding to a commensurate 
solid film. In addition, one finds two coverage ranges, 
separated by a cusp, where Oe{0)^ shown in Figure O has 
a positive second derivative (meaning that these config- 
urations satisfy the condition of positive compressibility, 
djj./dO=(P[Oe{0)]/dO'^ > 0, and are thus thermodynam- 
ically stable); we will call the first of these IC-I, and 
the second IC-II. IC-I extends from l9=0.07598 A-^ to 
6>=0.08030 A-2, while IC-II extends from i9=0.08030 A-^ 
to ^=0.8635 A~^. The much richer behavior of e{0) in 
the case of 0-D2, is a direct consequence of the greater 
mass, and reduced zero-point motion, of this isotope; as 
a result, the potential energy plays a greater role in shap- 
ing the physical properties. Phases of intermediate cov- 
erages found for 0-D2 are apparently "washed out" by 
zero-point motion in the case of P-H2. While this re- 
sult is surprising, similar anomalous behavior is noted in 
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FIG. 3: (color online) The product Oe{0) (in KA-^) for p-Hs 
(o), and for 0-D2 (•) computed by PIGS (T=0), as a function 
of the coverage (in A~^). ▲ denote coverages within regions 
of positive curvature. 
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FIG. 4: (color online) Density profile of 0-D2 adsorbed to 
the composite substrate for 6'=0.08389 A"^ and 6>=0.07253 
A~^ at T=l K (solid line and dashed line respectively). The 
square represents the position of graphite (at z=0 A) and the 
circle represents the position of the Kr layer (at Z—3A6 A). 



related studies i i 

Of particular interest for the present study is IC-II, 
where the 0-D2 is found to form a structure best char- 
acterized as being composed of several crystalline do- 
mains commensurately occupying one of the two trian- 
gular sublattices available. Between these commensurate 
domains are regions of high 0-D2 density - domain walls. 
Molecules that are part of these domain walls have a de- 
creased average interparticle distance. One may expect 
that such 0-D2, having the same zero-point motion as 0- 



D2 in the commensurate domains, would have a greater 
likelihood of quantum exchange, in particular if there is 
domain-wall disorder. Indeed, IC-II includes precisely 
the coverages reported in the experiments of Ref. [HI to 
exhibit anomalous fluid-like signatures down to T=1.5 K, 
well below the expected freezing point. 

In order to provide accurate permutation statistics and 
superfluid density estimates, additional simulations were 
carried out for several coverages at finite temperature 
using the worm algorithm. Figure [21 inset, shows the 
energy values obtained in the range of temperature 0.25 
K < T < 5 K, for i9=0.07253 A-^ and i9=0.08389 A-^. 

Shown in Figure [H are vertical density profiles, p{z), 
for 19=0.07253 A-^ and (9=0.08389 A-^ at T=l K. One 
can see that the vertical extent of the 0-D2 film is slightly 
enhanced for increased coverages. 

While dynamical information cannot be extracted us- 
ing the computational techniques employed in this study, 
one does have access to real-space instantaneous config- 
urations of the system, from which one can extract qual- 
itative insight, and potentially use to preclude particu- 
lar system behaviors. Such representative many-particle 
configurations are shown in Figure [5] for 0-D2 and p- 
H2 with 19=0.07253 A-^ at T = 0, and for 0-D2 with 
19=0.08389 A-2 at T=0.5 K. In each of (a) to (c), we 
find a high degree of particle localization, even along 
the domain walls observed in (c). No single configura- 
tion observed was found to be inconsistent with the in- 
stantaneous configurations one would necessarily observe 
in a domain-wall fluid, while the spatial configuration 
of domain-walls was not fixed as the system evolved in 
Monte Carlo time. 

Finally, consistent with the high degree of 0-D2 lo- 
calization observed, permutations are found to not take 
place at any temperature for any coverage examined, and, 
accordingly, the finite temperature estimator for the su- 
perfluid density consistently yields a value of Ps=0. In 
light of the absence of structural and energetic differences 
found between T=0 and T=0.25 K, it is unlikely that 
novel phases of 0-D2 (namely, phases with a finite super- 
fluid density) remain unexplored at intermediate temper- 
atures. 



V. DISCUSSION 

Using numerically exact ground state Quantum Monte 
Carlo (PIGS) and finite-temperature path integral Monte 
Carlo (worm algorithm) methods, we studied ^H2 and 
0-D2 adsorption onto a graphite substrate preplated with 
a compressed krypton monolayer. We performed calcula- 
tions based a simple model, in which graphite corrugation 
is ignored, the Kr atoms in the spacer layer are assumed 
static and point-like, and p-H2(o-D2)-substrate interac- 
tions are given by Lennard- Jones type potentials. Using 
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FIG. 5: (color online) Snapshots of typical configurations of 
P-II2 molecules adsorbed to the graphite/compressed Kr sub- 
strate. The positions of all molecules at each one of the imag- 
inary time slices are shown as discrete paths. Distances are 
expressed in A. These are low resolution images; contact the 
authors for high resolution images. 



an exact ground state technique, for both p-B.2 and 0- 
D2, we find that a solid monolayer commensurate with 
the Kr layer is thermodynamically stable, and, in both 
cases, is the coverage corresponding to the minimum en- 
ergy per particle. For 0-D2, we also find that there are 
two distinct compressible incommensurate solid regions, 
with the denser region, IC-II, corresponding to regions 
of commensuration separated by domain walls. Investi- 
gation of this coverage region at finite temperature, us- 
ing the worm algorithm, down to T=0.25 K, yields no 
evidence of superfluidity, consistent with our observa- 
tion of a high degree of localization of 0-D2 molecules. 
No quantum exchanges, and by extension permutations, 
are found to take place, though they are sampled effi- 
ciently. Though our simulations support the notion of a 
"domain- wall fluid", as all instantaneous configurations 
generated by our QMC algorithm do indeed show such 
domain walls, this mechanism is found not to support 
superfluidity of either hydrogen isotope. 

There are obviously several sources of uncertainty in 
this calculation which need to be discussed. The poten- 
tials used to describe the interactions between the P-H2 
and the substrate are very rough; this does not seem 
too important an issue, as far as the interaction of P-H2 
molecules with graphite is concerned, given the relatively 
large average distance at which molecules sit, supported 
by calculations using explicit corrugation of the graphite 
substrate. On the other hand, a more realistic interaction 
potential between P-H2 and krypton may quantitatively 
alter the energetics shown here. Despite these issues, 
and other simplifications, it does not seem likely that the 
qualitative structural information will change. Thus, we 
conclude this system is not a good candidate for further 
consideration in the search for superfluid 0-D2 or p-B.2' 

In general, the results shown here confirm those ob- 
tained in a other studies, to the effect that frustration 
arising from incommensuration does not lead to the for- 
mation of a stable superfluid phase of molecular hydrogen 
at low T. A more complex physical mechanism, causing 
a weakening of the effective interaction between two hy- 
drogen molecules may have to be investigated, if some 
progress is to be achieved in this direction. Possibly, 
the interaction of the hydrogen molecules with the sur- 
face electrons of a metal substrate, or of a nanostructure, 
might have such a "renormalizing" effect. 
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